/*This file is part of the FEBio source code and is licensed under the MIT license
listed below.

See Copyright-FEBio.txt for details.

Copyright (c) 2020 University of Utah, The Trustees of Columbia University in 
the City of New York, and others.

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.*/
#include "stdafx.h"
#include "FESolidElementShape.h"

//=============================================================================
//                                 H E X 8
//=============================================================================

//-----------------------------------------------------------------------------
void FEHex8::shape_fnc(double* H, double r, double s, double t)
{
	H[0] = 0.125*(1 - r)*(1 - s)*(1 - t);
	H[1] = 0.125*(1 + r)*(1 - s)*(1 - t);
	H[2] = 0.125*(1 + r)*(1 + s)*(1 - t);
	H[3] = 0.125*(1 - r)*(1 + s)*(1 - t);
	H[4] = 0.125*(1 - r)*(1 - s)*(1 + t);
	H[5] = 0.125*(1 + r)*(1 - s)*(1 + t);
	H[6] = 0.125*(1 + r)*(1 + s)*(1 + t);
	H[7] = 0.125*(1 - r)*(1 + s)*(1 + t);
}

//-----------------------------------------------------------------------------
//! values of shape function derivatives
void FEHex8::shape_deriv(double* Hr, double* Hs, double* Ht, double r, double s, double t)
{
	Hr[0] = -0.125*(1 - s)*(1 - t);
	Hr[1] = 0.125*(1 - s)*(1 - t);
	Hr[2] = 0.125*(1 + s)*(1 - t);
	Hr[3] = -0.125*(1 + s)*(1 - t);
	Hr[4] = -0.125*(1 - s)*(1 + t);
	Hr[5] = 0.125*(1 - s)*(1 + t);
	Hr[6] = 0.125*(1 + s)*(1 + t);
	Hr[7] = -0.125*(1 + s)*(1 + t);

	Hs[0] = -0.125*(1 - r)*(1 - t);
	Hs[1] = -0.125*(1 + r)*(1 - t);
	Hs[2] = 0.125*(1 + r)*(1 - t);
	Hs[3] = 0.125*(1 - r)*(1 - t);
	Hs[4] = -0.125*(1 - r)*(1 + t);
	Hs[5] = -0.125*(1 + r)*(1 + t);
	Hs[6] = 0.125*(1 + r)*(1 + t);
	Hs[7] = 0.125*(1 - r)*(1 + t);

	Ht[0] = -0.125*(1 - r)*(1 - s);
	Ht[1] = -0.125*(1 + r)*(1 - s);
	Ht[2] = -0.125*(1 + r)*(1 + s);
	Ht[3] = -0.125*(1 - r)*(1 + s);
	Ht[4] = 0.125*(1 - r)*(1 - s);
	Ht[5] = 0.125*(1 + r)*(1 - s);
	Ht[6] = 0.125*(1 + r)*(1 + s);
	Ht[7] = 0.125*(1 - r)*(1 + s);
}

//-----------------------------------------------------------------------------
//! values of shape function second derivatives
void FEHex8::shape_deriv2(double* Hrr, double* Hss, double* Htt, double* Hrs, double* Hst, double* Hrt, double r, double s, double t)
{
	Hrr[0] = 0.0; Hss[0] = 0.0; Htt[0] = 0.0;
	Hrr[1] = 0.0; Hss[1] = 0.0; Htt[1] = 0.0;
	Hrr[2] = 0.0; Hss[2] = 0.0; Htt[2] = 0.0;
	Hrr[3] = 0.0; Hss[3] = 0.0; Htt[3] = 0.0;
	Hrr[4] = 0.0; Hss[4] = 0.0; Htt[4] = 0.0;
	Hrr[5] = 0.0; Hss[5] = 0.0; Htt[5] = 0.0;
	Hrr[6] = 0.0; Hss[6] = 0.0; Htt[6] = 0.0;
	Hrr[7] = 0.0; Hss[7] = 0.0; Htt[7] = 0.0;

	Hrs[0] = 0.125*(1 - t);
	Hrs[1] = -0.125*(1 - t);
	Hrs[2] = 0.125*(1 - t);
	Hrs[3] = -0.125*(1 - t);
	Hrs[4] = 0.125*(1 + t);
	Hrs[5] = -0.125*(1 + t);
	Hrs[6] = 0.125*(1 + t);
	Hrs[7] = -0.125*(1 + t);

	Hrt[0] = 0.125*(1 - s);
	Hrt[1] = -0.125*(1 - s);
	Hrt[2] = -0.125*(1 + s);
	Hrt[3] = 0.125*(1 + s);
	Hrt[4] = -0.125*(1 - s);
	Hrt[5] = 0.125*(1 - s);
	Hrt[6] = 0.125*(1 + s);
	Hrt[7] = -0.125*(1 + s);

	Hst[0] = 0.125*(1 - r);
	Hst[1] = 0.125*(1 + r);
	Hst[2] = -0.125*(1 + r);
	Hst[3] = -0.125*(1 - r);
	Hst[4] = -0.125*(1 - r);
	Hst[5] = -0.125*(1 + r);
	Hst[6] = 0.125*(1 + r);
	Hst[7] = 0.125*(1 - r);
}

//=============================================================================
//                                 T E T 4
//=============================================================================

//-----------------------------------------------------------------------------
//! values of shape functions
void FETet4::shape_fnc(double* H, double r, double s, double t)
{
	H[0] = 1 - r - s - t;
	H[1] = r;
	H[2] = s;
	H[3] = t;
}

//-----------------------------------------------------------------------------
//! values of shape function derivatives
void FETet4::shape_deriv(double* Hr, double* Hs, double* Ht, double r, double s, double t)
{
	Hr[0] = -1; Hs[0] = -1; Ht[0] = -1;
	Hr[1] = 1;	Hs[1] = 0; Ht[1] = 0;
	Hr[2] = 0;	Hs[2] = 1; Ht[2] = 0;
	Hr[3] = 0;	Hs[3] = 0; Ht[3] = 1;
}

//-----------------------------------------------------------------------------
//! values of shape function second derivatives
void FETet4::shape_deriv2(double* Hrr, double* Hss, double* Htt, double* Hrs, double* Hst, double* Hrt, double r, double s, double t)
{
	Hrr[0] = 0.0; Hss[0] = 0.0; Htt[0] = 0.0;
	Hrr[1] = 0.0; Hss[1] = 0.0; Htt[1] = 0.0;
	Hrr[2] = 0.0; Hss[2] = 0.0; Htt[2] = 0.0;
	Hrr[3] = 0.0; Hss[3] = 0.0; Htt[3] = 0.0;

	Hrs[0] = 0.0; Hst[0] = 0.0; Hrt[0] = 0.0;
	Hrs[1] = 0.0; Hst[1] = 0.0; Hrt[1] = 0.0;
	Hrs[2] = 0.0; Hst[2] = 0.0; Hrt[2] = 0.0;
	Hrs[3] = 0.0; Hst[3] = 0.0; Hrt[3] = 0.0;
}

//=============================================================================
//                                 T E T 5
//=============================================================================

//-----------------------------------------------------------------------------
//! values of shape functions
void FETet5::shape_fnc(double* H, double r, double s, double t)
{
	H[0] = 1 - r - s - t;
	H[1] = r;
	H[2] = s;
	H[3] = t;
	H[4] = 256.0 * H[0] * H[1] * H[2] * H[3];

	H[0] -= 0.25*H[4];
	H[1] -= 0.25*H[4];
	H[2] -= 0.25*H[4];
	H[3] -= 0.25*H[4];
}

//-----------------------------------------------------------------------------
//! values of shape function derivatives
void FETet5::shape_deriv(double* Hr, double* Hs, double* Ht, double r, double s, double t)
{
	Hr[0] = -1; Hs[0] = -1; Ht[0] = -1;
	Hr[1] = 1;	Hs[1] = 0; Ht[1] = 0;
	Hr[2] = 0;	Hs[2] = 1; Ht[2] = 0;
	Hr[3] = 0;	Hs[3] = 0; Ht[3] = 1;

	Hr[4] = 256.0*s*t*(1.0 - 2.0*r - s - t);
	Hs[4] = 256.0*r*t*(1.0 - r - 2.0*s - t);
	Ht[4] = 256.0*r*s*(1.0 - r - s - 2.0*t);

	Hr[0] -= 0.25*Hr[4]; Hr[1] -= 0.25*Hr[4]; Hr[2] -= 0.25*Hr[4]; Hr[3] -= 0.25*Hr[4];
	Hs[0] -= 0.25*Hs[4]; Hs[1] -= 0.25*Hs[4]; Hs[2] -= 0.25*Hs[4]; Hs[3] -= 0.25*Hs[4];
	Ht[0] -= 0.25*Ht[4]; Ht[1] -= 0.25*Ht[4]; Ht[2] -= 0.25*Ht[4]; Ht[3] -= 0.25*Ht[4];
}

//-----------------------------------------------------------------------------
//! values of shape function second derivatives
void FETet5::shape_deriv2(double* Hrr, double* Hss, double* Htt, double* Hrs, double* Hst, double* Hrt, double r, double s, double t)
{
	Hrr[4] = -512 * s*t;
	Hss[4] = -512 * r*t;
	Htt[4] = -512 * r*s;

	Hrs[4] = 256.0*t*(1.0 - 2.0*r - 2.0*s - t);
	Hst[4] = 256.0*r*(1.0 - r - 2.0*s - 2.0*t);
	Hrt[4] = 256.0*s*(1.0 - 2.0*r - s - 2.0*t);

	Hrr[0] = -0.25*Hrr[4]; Hss[0] = -0.25*Hss[4]; Htt[0] = -0.25*Htt[4];
	Hrr[1] = -0.25*Hrr[4]; Hss[1] = -0.25*Hss[4]; Htt[1] = -0.25*Htt[4];
	Hrr[2] = -0.25*Hrr[4]; Hss[2] = -0.25*Hss[4]; Htt[2] = -0.25*Htt[4];
	Hrr[3] = -0.25*Hrr[4]; Hss[3] = -0.25*Hss[4]; Htt[3] = -0.25*Htt[4];

	Hrs[0] = -0.25*Hrs[4]; Hst[0] = -0.25*Hst[4]; Hrt[0] = -Hrt[4];
	Hrs[1] = -0.25*Hrs[4]; Hst[1] = -0.25*Hst[4]; Hrt[1] = -Hrt[4];
	Hrs[2] = -0.25*Hrs[4]; Hst[2] = -0.25*Hst[4]; Hrt[2] = -Hrt[4];
	Hrs[3] = -0.25*Hrs[4]; Hst[3] = -0.25*Hst[4]; Hrt[3] = -Hrt[4];
}

//=============================================================================
//                                 P E N T A 6
//=============================================================================

//-----------------------------------------------------------------------------
//! values of shape functions
void FEPenta6::shape_fnc(double* H, double r, double s, double t)
{
	H[0] = 0.5*(1 - t)*(1 - r - s);
	H[1] = 0.5*(1 - t)*r;
	H[2] = 0.5*(1 - t)*s;
	H[3] = 0.5*(1 + t)*(1 - r - s);
	H[4] = 0.5*(1 + t)*r;
	H[5] = 0.5*(1 + t)*s;
}

//-----------------------------------------------------------------------------
//! values of shape function derivatives
void FEPenta6::shape_deriv(double* Hr, double* Hs, double* Ht, double r, double s, double t)
{
	Hr[0] = -0.5*(1 - t);
	Hr[1] = 0.5*(1 - t);
	Hr[2] = 0.0;
	Hr[3] = -0.5*(1 + t);
	Hr[4] = 0.5*(1 + t);
	Hr[5] = 0.0;

	Hs[0] = -0.5*(1 - t);
	Hs[1] = 0.0;
	Hs[2] = 0.5*(1 - t);
	Hs[3] = -0.5*(1 + t);
	Hs[4] = 0.0;
	Hs[5] = 0.5*(1 + t);

	Ht[0] = -0.5*(1 - r - s);
	Ht[1] = -0.5*r;
	Ht[2] = -0.5*s;
	Ht[3] = 0.5*(1 - r - s);
	Ht[4] = 0.5*r;
	Ht[5] = 0.5*s;
}

//-----------------------------------------------------------------------------
//! values of shape function second derivatives
void FEPenta6::shape_deriv2(double* Hrr, double* Hss, double* Htt, double* Hrs, double* Hst, double* Hrt, double r, double s, double t)
{
	Hrr[0] = 0.0; Hss[0] = 0.0; Htt[0] = 0.0;
	Hrr[1] = 0.0; Hss[1] = 0.0; Htt[1] = 0.0;
	Hrr[2] = 0.0; Hss[2] = 0.0; Htt[2] = 0.0;
	Hrr[3] = 0.0; Hss[3] = 0.0; Htt[3] = 0.0;
	Hrr[4] = 0.0; Hss[4] = 0.0; Htt[4] = 0.0;
	Hrr[5] = 0.0; Hss[5] = 0.0; Htt[5] = 0.0;

	Hrs[0] = 0.0; Hst[0] = 0.5; Hrt[0] = 0.5;
	Hrs[1] = 0.0; Hst[1] = 0.0; Hrt[1] = -0.5;
	Hrs[2] = 0.0; Hst[2] = -0.5; Hrt[2] = 0.0;
	Hrs[3] = 0.0; Hst[3] = -0.5; Hrt[3] = -0.5;
	Hrs[4] = 0.0; Hst[4] = 0.0; Hrt[4] = 0.5;
	Hrs[5] = 0.0; Hst[5] = 0.5; Hrt[5] = 0.0;
}

//=============================================================================
//                                 P E N T A 1 5
//=============================================================================

//-----------------------------------------------------------------------------
//! values of shape functions
void FEPenta15::shape_fnc(double* H, double r, double s, double t)
{
	double u = 1 - r - s;

	H[0] = -((1 - t*t)*u) / 2. + ((1 - t)*u*(-1 + 2 * u)) / 2.;
	H[1] = (r*(-1 + 2 * r)*(1 - t)) / 2. - (r*(1 - t*t)) / 2.;
	H[2] = (s*(-1 + 2 * s)*(1 - t)) / 2. - (s*(1 - t*t)) / 2.;
	H[3] = -((1 - t*t)*u) / 2. + ((1 + t)*u*(-1 + 2 * u)) / 2.;
	H[4] = (r*(-1 + 2 * r)*(1 + t)) / 2. - (r*(1 - t*t)) / 2.;
	H[5] = (s*(-1 + 2 * s)*(1 + t)) / 2. - (s*(1 - t*t)) / 2.;
	H[6] = 2 * r*(1 - t)*u;
	H[7] = 2 * r*s*(1 - t);
	H[8] = 2 * s*(1 - t)*u;
	H[9] = 2 * r*(1 + t)*u;
	H[10] = 2 * r*s*(1 + t);
	H[11] = 2 * s*(1 + t)*u;
	H[12] = (1 - t*t)*u;
	H[13] = r*(1 - t*t);
	H[14] = s*(1 - t*t);
}

//-----------------------------------------------------------------------------
//! values of shape function derivatives
void FEPenta15::shape_deriv(double* Hr, double* Hs, double* Ht, double r, double s, double t)
{
	Hr[0] = -((-1 + t)*(-2 + 4 * r + 4 * s + t)) / 2.;
	Hr[1] = (-2 - 4 * r*(-1 + t) + t + t*t) / 2.;
	Hr[2] = 0;
	Hr[3] = ((-2 + 4 * r + 4 * s - t)*(1 + t)) / 2.;
	Hr[4] = ((1 + t)*(-2 + 4 * r + t)) / 2.;
	Hr[5] = 0;
	Hr[6] = 2 * (-1 + 2 * r + s)*(-1 + t);
	Hr[7] = -2 * s*(-1 + t);
	Hr[8] = 2 * s*(-1 + t);
	Hr[9] = -2 * (-1 + 2 * r + s)*(1 + t);
	Hr[10] = 2 * s*(1 + t);
	Hr[11] = -2 * s*(1 + t);
	Hr[12] = -1 + t*t;
	Hr[13] = 1 - t*t;
	Hr[14] = 0;

	Hs[0] = -((-1 + t)*(-2 + 4 * r + 4 * s + t)) / 2.;
	Hs[1] = 0;
	Hs[2] = (-2 - 4 * s*(-1 + t) + t + t*t) / 2.;
	Hs[3] = ((-2 + 4 * r + 4 * s - t)*(1 + t)) / 2.;
	Hs[4] = 0;
	Hs[5] = ((1 + t)*(-2 + 4 * s + t)) / 2.;
	Hs[6] = 2 * r*(-1 + t);
	Hs[7] = -2 * r*(-1 + t);
	Hs[8] = 2 * (-1 + r + 2 * s)*(-1 + t);
	Hs[9] = -2 * r*(1 + t);
	Hs[10] = 2 * r*(1 + t);
	Hs[11] = -2 * (-1 + r + 2 * s)*(1 + t);
	Hs[12] = -1 + t*t;
	Hs[13] = 0;
	Hs[14] = 1 - t*t;

	Ht[0] = -((-1 + r + s)*(-1 + 2 * r + 2 * s + 2 * t)) / 2.;
	Ht[1] = (r*(1 - 2 * r + 2 * t)) / 2.;
	Ht[2] = (s*(1 - 2 * s + 2 * t)) / 2.;
	Ht[3] = ((-1 + r + s)*(-1 + 2 * r + 2 * s - 2 * t)) / 2.;
	Ht[4] = r*(-0.5 + r + t);
	Ht[5] = s*(-0.5 + s + t);
	Ht[6] = 2 * r*(-1 + r + s);
	Ht[7] = -2 * r*s;
	Ht[8] = 2 * s*(-1 + r + s);
	Ht[9] = -2 * r*(-1 + r + s);
	Ht[10] = 2 * r*s;
	Ht[11] = -2 * s*(-1 + r + s);
	Ht[12] = 2 * (-1 + r + s)*t;
	Ht[13] = -2 * r*t;
	Ht[14] = -2 * s*t;
}

//-----------------------------------------------------------------------------
//! values of shape function second derivatives
void FEPenta15::shape_deriv2(double* Hrr, double* Hss, double* Htt, double* Hrs, double* Hst, double* Hrt, double r, double s, double t)
{
	Hrr[0] = 2 - 2 * t;
	Hrr[1] = 2 - 2 * t;
	Hrr[2] = 0;
	Hrr[3] = 2 * (1 + t);
	Hrr[4] = 2 * (1 + t);
	Hrr[5] = 0;
	Hrr[6] = 4 * (-1 + t);
	Hrr[7] = 0;
	Hrr[8] = 0;
	Hrr[9] = -4 * (1 + t);
	Hrr[10] = 0;
	Hrr[11] = 0;
	Hrr[12] = 0;
	Hrr[13] = 0;
	Hrr[14] = 0;

	Hss[0] = 2 - 2 * t;
	Hss[1] = 0;
	Hss[2] = 2 - 2 * t;
	Hss[3] = 2 * (1 + t);
	Hss[4] = 0;
	Hss[5] = 2 * (1 + t);
	Hss[6] = 0;
	Hss[7] = 0;
	Hss[8] = 4 * (-1 + t);
	Hss[9] = 0;
	Hss[10] = 0;
	Hss[11] = -4 * (1 + t);
	Hss[12] = 0;
	Hss[13] = 0;
	Hss[14] = 0;

	Htt[0] = 1 - r - s;
	Htt[1] = r;
	Htt[2] = s;
	Htt[3] = 1 - r - s;
	Htt[4] = r;
	Htt[5] = s;
	Htt[6] = 0;
	Htt[7] = 0;
	Htt[8] = 0;
	Htt[9] = 0;
	Htt[10] = 0;
	Htt[11] = 0;
	Htt[12] = 2 * (-1 + r + s);
	Htt[13] = -2 * r;
	Htt[14] = -2 * s;

	Hrs[0] = 2 - 2 * t;
	Hrs[1] = 0;
	Hrs[2] = 0;
	Hrs[3] = 2 * (1 + t);
	Hrs[4] = 0;
	Hrs[5] = 0;
	Hrs[6] = 2 * (-1 + t);
	Hrs[7] = 2 - 2 * t;
	Hrs[8] = 2 * (-1 + t);
	Hrs[9] = -2 * (1 + t);
	Hrs[10] = 2 * (1 + t);
	Hrs[11] = -2 * (1 + t);
	Hrs[12] = 0;
	Hrs[13] = 0;
	Hrs[14] = 0;

	Hst[0] = 1.5 - 2 * r - 2 * s - t;
	Hst[1] = 0;
	Hst[2] = 0.5 - 2 * s + t;
	Hst[3] = -1.5 + 2 * r + 2 * s - t;
	Hst[4] = 0;
	Hst[5] = -0.5 + 2 * s + t;
	Hst[6] = 2 * r;
	Hst[7] = -2 * r;
	Hst[8] = 2 * (-1 + r + 2 * s);
	Hst[9] = -2 * r;
	Hst[10] = 2 * r;
	Hst[11] = -2 * (-1 + r + 2 * s);
	Hst[12] = 2 * t;
	Hst[13] = 0;
	Hst[14] = -2 * t;

	Hrt[0] = 1.5 - 2 * r - 2 * s - t;
	Hrt[1] = 0.5 - 2 * r + t;
	Hrt[2] = 0;
	Hrt[3] = -1.5 + 2 * r + 2 * s - t;
	Hrt[4] = -0.5 + 2 * r + t;
	Hrt[5] = 0;
	Hrt[6] = 2 * (-1 + 2 * r + s);
	Hrt[7] = -2 * s;
	Hrt[8] = 2 * s;
	Hrt[9] = -2 * (-1 + 2 * r + s);
	Hrt[10] = 2 * s;
	Hrt[11] = -2 * s;
	Hrt[12] = 2 * t;
	Hrt[13] = -2 * t;
	Hrt[14] = 0;
}

//=============================================================================
//                                 T E T 1 0
//=============================================================================

//-----------------------------------------------------------------------------
//! values of shape functions
void FETet10::shape_fnc(double* H, double r, double s, double t)
{
	double r1 = 1.0 - r - s - t;
	double r2 = r;
	double r3 = s;
	double r4 = t;

	H[0] = r1*(2.0*r1 - 1.0);
	H[1] = r2*(2.0*r2 - 1.0);
	H[2] = r3*(2.0*r3 - 1.0);
	H[3] = r4*(2.0*r4 - 1.0);
	H[4] = 4.0*r1*r2;
	H[5] = 4.0*r2*r3;
	H[6] = 4.0*r3*r1;
	H[7] = 4.0*r1*r4;
	H[8] = 4.0*r2*r4;
	H[9] = 4.0*r3*r4;
}

//-----------------------------------------------------------------------------
//! values of shape function derivatives
void FETet10::shape_deriv(double* Hr, double* Hs, double* Ht, double r, double s, double t)
{
	Hr[0] = -3.0 + 4.0*r + 4.0*(s + t);
	Hr[1] = 4.0*r - 1.0;
	Hr[2] = 0.0;
	Hr[3] = 0.0;
	Hr[4] = 4.0 - 8.0*r - 4.0*(s + t);
	Hr[5] = 4.0*s;
	Hr[6] = -4.0*s;
	Hr[7] = -4.0*t;
	Hr[8] = 4.0*t;
	Hr[9] = 0.0;

	Hs[0] = -3.0 + 4.0*s + 4.0*(r + t);
	Hs[1] = 0.0;
	Hs[2] = 4.0*s - 1.0;
	Hs[3] = 0.0;
	Hs[4] = -4.0*r;
	Hs[5] = 4.0*r;
	Hs[6] = 4.0 - 8.0*s - 4.0*(r + t);
	Hs[7] = -4.0*t;
	Hs[8] = 0.0;
	Hs[9] = 4.0*t;

	Ht[0] = -3.0 + 4.0*t + 4.0*(r + s);
	Ht[1] = 0.0;
	Ht[2] = 0.0;
	Ht[3] = 4.0*t - 1.0;
	Ht[4] = -4.0*r;
	Ht[5] = 0.0;
	Ht[6] = -4.0*s;
	Ht[7] = 4.0 - 8.0*t - 4.0*(r + s);
	Ht[8] = 4.0*r;
	Ht[9] = 4.0*s;
}

//-----------------------------------------------------------------------------
//! values of shape function second derivatives
void FETet10::shape_deriv2(double* Hrr, double* Hss, double* Htt, double* Hrs, double* Hst, double* Hrt, double r, double s, double t)
{
	Hrr[0] = 4.0; Hss[0] = 4.0; Htt[0] = 4.0;
	Hrr[1] = 4.0; Hss[1] = 0.0; Htt[1] = 0.0;
	Hrr[2] = 0.0; Hss[2] = 4.0; Htt[2] = 0.0;
	Hrr[3] = 0.0; Hss[3] = 0.0; Htt[3] = 4.0;
	Hrr[4] = -8.0; Hss[4] = 0.0; Htt[4] = 0.0;
	Hrr[5] = 0.0; Hss[5] = 0.0; Htt[5] = 0.0;
	Hrr[6] = 0.0; Hss[6] = -8.0; Htt[6] = 0.0;
	Hrr[7] = 0.0; Hss[7] = 0.0; Htt[7] = -8.0;
	Hrr[8] = 0.0; Hss[8] = 0.0; Htt[8] = 0.0;
	Hrr[9] = 0.0; Hss[9] = 0.0; Htt[9] = 0.0;

	Hrs[0] = 4.0; Hst[0] = 4.0; Hrt[0] = 4.0;
	Hrs[1] = 0.0; Hst[1] = 0.0; Hrt[1] = 0.0;
	Hrs[2] = 0.0; Hst[2] = 0.0; Hrt[2] = 0.0;
	Hrs[3] = 0.0; Hst[3] = 0.0; Hrt[3] = 0.0;
	Hrs[4] = -4.0; Hst[4] = 0.0; Hrt[4] = -4.0;
	Hrs[5] = 4.0; Hst[5] = 0.0; Hrt[5] = 0.0;
	Hrs[6] = -4.0; Hst[6] = -4.0; Hrt[6] = 0.0;
	Hrs[7] = 0.0; Hst[7] = -4.0; Hrt[7] = -4.0;
	Hrs[8] = 0.0; Hst[8] = 0.0; Hrt[8] = 4.0;
	Hrs[9] = 0.0; Hst[9] = 4.0; Hrt[9] = 0.0;
}


//=============================================================================
//                           T E T 1 5
//=============================================================================

//-----------------------------------------------------------------------------
void FETet15::shape_fnc(double* H, double r, double s, double t)
{
	double r1 = 1.0 - r - s - t;
	double r2 = r;
	double r3 = s;
	double r4 = t;

	H[14] = 256 * r1*r2*r3*r4;

	H[10] = 27.0*r1*r2*r3 - 27.0*H[14] / 64.0;
	H[11] = 27.0*r1*r2*r4 - 27.0*H[14] / 64.0;
	H[12] = 27.0*r2*r3*r4 - 27.0*H[14] / 64.0;
	H[13] = 27.0*r3*r1*r4 - 27.0*H[14] / 64.0;

	H[0] = r1*(2.0*r1 - 1.0) + (H[10] + H[11] + H[13]) / 9.0 + H[14] / 8.0;
	H[1] = r2*(2.0*r2 - 1.0) + (H[10] + H[11] + H[12]) / 9.0 + H[14] / 8.0;
	H[2] = r3*(2.0*r3 - 1.0) + (H[10] + H[12] + H[13]) / 9.0 + H[14] / 8.0;
	H[3] = r4*(2.0*r4 - 1.0) + (H[11] + H[12] + H[13]) / 9.0 + H[14] / 8.0;

	H[4] = 4.0*r1*r2 - 4.0*(H[10] + H[11]) / 9.0 - H[14] / 4.0;
	H[5] = 4.0*r2*r3 - 4.0*(H[10] + H[12]) / 9.0 - H[14] / 4.0;
	H[6] = 4.0*r3*r1 - 4.0*(H[10] + H[13]) / 9.0 - H[14] / 4.0;
	H[7] = 4.0*r1*r4 - 4.0*(H[11] + H[13]) / 9.0 - H[14] / 4.0;
	H[8] = 4.0*r2*r4 - 4.0*(H[11] + H[12]) / 9.0 - H[14] / 4.0;
	H[9] = 4.0*r3*r4 - 4.0*(H[12] + H[13]) / 9.0 - H[14] / 4.0;
}

//-----------------------------------------------------------------------------
void FETet15::shape_deriv(double* Hr, double* Hs, double* Ht, double r, double s, double t)
{
	double u = 1.0 - r - s - t;

	Hr[14] = 256.0*s*t*(u - r);
	Hs[14] = 256.0*r*t*(u - s);
	Ht[14] = 256.0*r*s*(u - t);

	Hr[10] = 27.0*s*(u - r) - 27.0*Hr[14] / 64.0;
	Hr[11] = 27.0*t*(u - r) - 27.0*Hr[14] / 64.0;
	Hr[12] = 27.0*s*t - 27.0*Hr[14] / 64.0;
	Hr[13] = -27.0*s*t - 27.0*Hr[14] / 64.0;

	Hs[10] = 27.0*r*(u - s) - 27.0*Hs[14] / 64.0;
	Hs[11] = -27.0*r*t - 27.0*Hs[14] / 64.0;
	Hs[12] = 27.0*r*t - 27.0*Hs[14] / 64.0;
	Hs[13] = 27.0*t*(u - s) - 27.0*Hs[14] / 64.0;

	Ht[10] = -27.0*r*s - 27.0*Ht[14] / 64.0;
	Ht[11] = 27.0*r*(u - t) - 27.0*Ht[14] / 64.0;
	Ht[12] = 27.0*r*s - 27.0*Ht[14] / 64.0;
	Ht[13] = 27.0*s*(u - t) - 27.0*Ht[14] / 64.0;

	Hr[0] = -(4.0*u - 1.0) + (Hr[10] + Hr[11] + Hr[13]) / 9.0 + Hr[14] / 8.0;
	Hr[1] = (4.0*r - 1.0) + (Hr[10] + Hr[11] + Hr[12]) / 9.0 + Hr[14] / 8.0;
	Hr[2] = 0.0 + (Hr[10] + Hr[12] + Hr[13]) / 9.0 + Hr[14] / 8.0;
	Hr[3] = 0.0 + (Hr[11] + Hr[12] + Hr[13]) / 9.0 + Hr[14] / 8.0;
	Hr[4] = 4.0*(u - r) - 4.0*(Hr[10] + Hr[11]) / 9.0 - Hr[14] / 4.0;
	Hr[5] = 4.0*s - 4.0*(Hr[10] + Hr[12]) / 9.0 - Hr[14] / 4.0;
	Hr[6] = -4.0*s - 4.0*(Hr[10] + Hr[13]) / 9.0 - Hr[14] / 4.0;
	Hr[7] = -4.0*t - 4.0*(Hr[11] + Hr[13]) / 9.0 - Hr[14] / 4.0;
	Hr[8] = 4.0*t - 4.0*(Hr[11] + Hr[12]) / 9.0 - Hr[14] / 4.0;
	Hr[9] = 0.0 - 4.0*(Hr[12] + Hr[13]) / 9.0 - Hr[14] / 4.0;

	Hs[0] = -(4.0*u - 1.0) + (Hs[10] + Hs[11] + Hs[13]) / 9.0 + Hs[14] / 8.0;
	Hs[1] = 0.0 + (Hs[10] + Hs[11] + Hs[12]) / 9.0 + Hs[14] / 8.0;
	Hs[2] = (4.0*s - 1.0) + (Hs[10] + Hs[12] + Hs[13]) / 9.0 + Hs[14] / 8.0;
	Hs[3] = 0.0 + (Hs[11] + Hs[12] + Hs[13]) / 9.0 + Hs[14] / 8.0;
	Hs[4] = -4.0*r - 4.0*(Hs[10] + Hs[11]) / 9.0 - Hs[14] / 4.0;
	Hs[5] = 4.0*r - 4.0*(Hs[10] + Hs[12]) / 9.0 - Hs[14] / 4.0;
	Hs[6] = 4.0*(u - s) - 4.0*(Hs[10] + Hs[13]) / 9.0 - Hs[14] / 4.0;
	Hs[7] = -4.0*t - 4.0*(Hs[11] + Hs[13]) / 9.0 - Hs[14] / 4.0;
	Hs[8] = 0.0 - 4.0*(Hs[11] + Hs[12]) / 9.0 - Hs[14] / 4.0;
	Hs[9] = 4.0*t - 4.0*(Hs[12] + Hs[13]) / 9.0 - Hs[14] / 4.0;

	Ht[0] = -(4.0*u - 1.0) + (Ht[10] + Ht[11] + Ht[13]) / 9.0 + Ht[14] / 8.0;
	Ht[1] = 0.0 + (Ht[10] + Ht[11] + Ht[12]) / 9.0 + Ht[14] / 8.0;
	Ht[2] = 0.0 + (Ht[10] + Ht[12] + Ht[13]) / 9.0 + Ht[14] / 8.0;
	Ht[3] = (4.0*t - 1.0) + (Ht[11] + Ht[12] + Ht[13]) / 9.0 + Ht[14] / 8.0;
	Ht[4] = -4.0*r - 4.0*(Ht[10] + Ht[11]) / 9.0 - Ht[14] / 4.0;
	Ht[5] = 0.0 - 4.0*(Ht[10] + Ht[12]) / 9.0 - Ht[14] / 4.0;
	Ht[6] = -4.0*s - 4.0*(Ht[10] + Ht[13]) / 9.0 - Ht[14] / 4.0;
	Ht[7] = 4.0*(u - t) - 4.0*(Ht[11] + Ht[13]) / 9.0 - Ht[14] / 4.0;
	Ht[8] = 4.0*r - 4.0*(Ht[11] + Ht[12]) / 9.0 - Ht[14] / 4.0;
	Ht[9] = 4.0*s - 4.0*(Ht[12] + Ht[13]) / 9.0 - Ht[14] / 4.0;
}

//-----------------------------------------------------------------------------
void FETet15::shape_deriv2(double* Hrr, double* Hss, double* Htt, double* Hrs, double* Hst, double* Hrt, double r, double s, double t)
{
	double u = 1.0 - r - s - t;

	Hrr[14] = -512 * s*t;
	Hss[14] = -512 * r*t;
	Htt[14] = -512 * r*s;
	Hrs[14] = 256.0*t*(u - r - s);
	Hst[14] = 256.0*r*(u - s - t);
	Hrt[14] = 256.0*s*(u - r - t);

	Hrr[10] = -54.0*s - 27.0*Hrr[14] / 64.0;
	Hss[10] = -54.0*r - 27.0*Hss[14] / 64.0;
	Htt[10] = 0.0 - 27.0*Htt[14] / 64.0;
	Hrs[10] = 27.0*(u - r - s) - 27.0*Hrs[14] / 64.0;
	Hst[10] = -27.0*r - 27.0*Hst[14] / 64.0;
	Hrt[10] = -27.0*s - 27.0*Hrt[14] / 64.0;

	Hrr[11] = -54.0*t - 27.0*Hrr[14] / 64.0;
	Hss[11] = 0.0 - 27.0*Hss[14] / 64.0;
	Htt[11] = -54.*r - 27.0*Htt[14] / 64.0;
	Hrs[11] = -27.0*t - 27.0*Hrs[14] / 64.0;
	Hst[11] = -27.0*r - 27.0*Hst[14] / 64.0;
	Hrt[11] = 27.0*(u - r - t) - 27.0*Hrt[14] / 64.0;

	Hrr[12] = 0.0 - 27.0*Hrr[14] / 64.0;
	Hss[12] = 0.0 - 27.0*Hss[14] / 64.0;
	Htt[12] = 0.0 - 27.0*Htt[14] / 64.0;
	Hrs[12] = 27.0*t - 27.0*Hrs[14] / 64.0;
	Hst[12] = 27.0*r - 27.0*Hst[14] / 64.0;
	Hrt[12] = 27.0*s - 27.0*Hrt[14] / 64.0;

	Hrr[13] = 0.0 - 27.0*Hrr[14] / 64.0;
	Hss[13] = -54.0*t - 27.0*Hss[14] / 64.0;
	Htt[13] = -54.0*s - 27.0*Htt[14] / 64.0;
	Hrs[13] = -27.0*t - 27.0*Hrs[14] / 64.0;
	Hst[13] = 27.0*(u - t - s) - 27.0*Hst[14] / 64.0;
	Hrt[13] = -27.0*s - 27.0*Hrt[14] / 64.0;

	Hrr[0] = 4.0 + (Hrr[10] + Hrr[11] + Hrr[13]) / 9.0 + Hrr[14] / 8.0;
	Hss[0] = 4.0 + (Hss[10] + Hss[11] + Hss[13]) / 9.0 + Hss[14] / 8.0;
	Htt[0] = 4.0 + (Htt[10] + Htt[11] + Htt[13]) / 9.0 + Htt[14] / 8.0;
	Hrs[0] = 4.0 + (Hrs[10] + Hrs[11] + Hrs[13]) / 9.0 + Hrs[14] / 8.0;
	Hst[0] = 4.0 + (Hst[10] + Hst[11] + Hst[13]) / 9.0 + Hst[14] / 8.0;
	Hrt[0] = 4.0 + (Hrt[10] + Hrt[11] + Hrt[13]) / 9.0 + Hrt[14] / 8.0;

	Hrr[1] = 4.0 + (Hrr[10] + Hrr[11] + Hrr[12]) / 9.0 + Hrr[14] / 8.0;
	Hss[1] = 0.0 + (Hss[10] + Hss[11] + Hss[12]) / 9.0 + Hss[14] / 8.0;
	Htt[1] = 0.0 + (Htt[10] + Htt[11] + Htt[12]) / 9.0 + Htt[14] / 8.0;
	Hrs[1] = 0.0 + (Hrs[10] + Hrs[11] + Hrs[12]) / 9.0 + Hrs[14] / 8.0;
	Hst[1] = 0.0 + (Hst[10] + Hst[11] + Hst[12]) / 9.0 + Hst[14] / 8.0;
	Hrt[1] = 0.0 + (Hrt[10] + Hrt[11] + Hrt[12]) / 9.0 + Hrt[14] / 8.0;

	Hrr[2] = 0.0 + (Hrr[10] + Hrr[12] + Hrr[13]) / 9.0 + Hrr[14] / 8.0;
	Hss[2] = 4.0 + (Hss[10] + Hss[12] + Hss[13]) / 9.0 + Hss[14] / 8.0;
	Htt[2] = 0.0 + (Htt[10] + Htt[12] + Htt[13]) / 9.0 + Htt[14] / 8.0;
	Hrs[2] = 0.0 + (Hrs[10] + Hrs[12] + Hrs[13]) / 9.0 + Hrs[14] / 8.0;
	Hst[2] = 0.0 + (Hst[10] + Hst[12] + Hst[13]) / 9.0 + Hst[14] / 8.0;
	Hrt[2] = 0.0 + (Hrt[10] + Hrt[12] + Hrt[13]) / 9.0 + Hrt[14] / 8.0;

	Hrr[3] = 0.0 + (Hrr[11] + Hrr[12] + Hrr[13]) / 9.0 + Hrr[14] / 8.0;
	Hss[3] = 0.0 + (Hss[11] + Hss[12] + Hss[13]) / 9.0 + Hss[14] / 8.0;
	Htt[3] = 4.0 + (Htt[11] + Htt[12] + Htt[13]) / 9.0 + Htt[14] / 8.0;
	Hrs[3] = 0.0 + (Hrs[11] + Hrs[12] + Hrs[13]) / 9.0 + Hrs[14] / 8.0;
	Hst[3] = 0.0 + (Hst[11] + Hst[12] + Hst[13]) / 9.0 + Hst[14] / 8.0;
	Hrt[3] = 0.0 + (Hrt[11] + Hrt[12] + Hrt[13]) / 9.0 + Hrt[14] / 8.0;

	Hrr[4] = -8.0 - 4.0*(Hrr[10] + Hrr[11]) / 9.0 - Hrr[14] / 4.0;
	Hss[4] = 0.0 - 4.0*(Hss[10] + Hss[11]) / 9.0 - Hss[14] / 4.0;
	Htt[4] = 0.0 - 4.0*(Htt[10] + Htt[11]) / 9.0 - Htt[14] / 4.0;
	Hrs[4] = -4.0 - 4.0*(Hrs[10] + Hrs[11]) / 9.0 - Hrs[14] / 4.0;
	Hst[4] = 0.0 - 4.0*(Hst[10] + Hst[11]) / 9.0 - Hst[14] / 4.0;
	Hrt[4] = -4.0 - 4.0*(Hrt[10] + Hrt[11]) / 9.0 - Hrt[14] / 4.0;

	Hrr[5] = 0.0 - 4.0*(Hrr[10] + Hrr[12]) / 9.0 - Hrr[14] / 4.0;
	Hss[5] = 0.0 - 4.0*(Hss[10] + Hss[12]) / 9.0 - Hss[14] / 4.0;
	Htt[5] = 0.0 - 4.0*(Htt[10] + Htt[12]) / 9.0 - Htt[14] / 4.0;
	Hrs[5] = 4.0 - 4.0*(Hrs[10] + Hrs[12]) / 9.0 - Hrs[14] / 4.0;
	Hst[5] = 0.0 - 4.0*(Hst[10] + Hst[12]) / 9.0 - Hst[14] / 4.0;
	Hrt[5] = 0.0 - 4.0*(Hrt[10] + Hrt[12]) / 9.0 - Hrt[14] / 4.0;

	Hrr[6] = 0.0 - 4.0*(Hrr[10] + Hrr[13]) / 9.0 - Hrr[14] / 4.0;
	Hss[6] = -8.0 - 4.0*(Hss[10] + Hss[13]) / 9.0 - Hss[14] / 4.0;
	Htt[6] = 0.0 - 4.0*(Htt[10] + Htt[13]) / 9.0 - Htt[14] / 4.0;
	Hrs[6] = -4.0 - 4.0*(Hrs[10] + Hrs[13]) / 9.0 - Hrs[14] / 4.0;
	Hst[6] = -4.0 - 4.0*(Hst[10] + Hst[13]) / 9.0 - Hst[14] / 4.0;
	Hrt[6] = 0.0 - 4.0*(Hrt[10] + Hrt[13]) / 9.0 - Hrt[14] / 4.0;

	Hrr[7] = 0.0 - 4.0*(Hrr[11] + Hrr[13]) / 9.0 - Hrr[14] / 4.0;
	Hss[7] = 0.0 - 4.0*(Hss[11] + Hss[13]) / 9.0 - Hss[14] / 4.0;
	Htt[7] = -8.0 - 4.0*(Htt[11] + Htt[13]) / 9.0 - Htt[14] / 4.0;
	Hrs[7] = 0.0 - 4.0*(Hrs[11] + Hrs[13]) / 9.0 - Hrs[14] / 4.0;
	Hst[7] = -4.0 - 4.0*(Hst[11] + Hst[13]) / 9.0 - Hst[14] / 4.0;
	Hrt[7] = -4.0 - 4.0*(Hrt[11] + Hrt[13]) / 9.0 - Hrt[14] / 4.0;

	Hrr[8] = 0.0 - 4.0*(Hrr[11] + Hrr[12]) / 9.0 - Hrr[14] / 4.0;
	Hss[8] = 0.0 - 4.0*(Hss[11] + Hss[12]) / 9.0 - Hss[14] / 4.0;
	Htt[8] = 0.0 - 4.0*(Htt[11] + Htt[12]) / 9.0 - Htt[14] / 4.0;
	Hrs[8] = 0.0 - 4.0*(Hrs[11] + Hrs[12]) / 9.0 - Hrs[14] / 4.0;
	Hst[8] = 0.0 - 4.0*(Hst[11] + Hst[12]) / 9.0 - Hst[14] / 4.0;
	Hrt[8] = 4.0 - 4.0*(Hrt[11] + Hrt[12]) / 9.0 - Hrt[14] / 4.0;

	Hrr[9] = 0.0 - 4.0*(Hrr[12] + Hrr[13]) / 9.0 - Hrr[14] / 4.0;
	Hss[9] = 0.0 - 4.0*(Hss[12] + Hss[13]) / 9.0 - Hss[14] / 4.0;
	Htt[9] = 0.0 - 4.0*(Htt[12] + Htt[13]) / 9.0 - Htt[14] / 4.0;
	Hrs[9] = 0.0 - 4.0*(Hrs[12] + Hrs[13]) / 9.0 - Hrs[14] / 4.0;
	Hst[9] = 4.0 - 4.0*(Hst[12] + Hst[13]) / 9.0 - Hst[14] / 4.0;
	Hrt[9] = 0.0 - 4.0*(Hrt[12] + Hrt[13]) / 9.0 - Hrt[14] / 4.0;
}

//=============================================================================
//                           T E T 2 0
//=============================================================================

//-----------------------------------------------------------------------------
//! values of shape functions
void FETet20::shape_fnc(double* H, double r, double s, double t)
{
	double L1 = 1.0 - r - s - t;
	double L2 = r;
	double L3 = s;
	double L4 = t;

	H[0] = 0.5*(3 * L1 - 1)*(3 * L1 - 2)*L1;
	H[1] = 0.5*(3 * L2 - 1)*(3 * L2 - 2)*L2;
	H[2] = 0.5*(3 * L3 - 1)*(3 * L3 - 2)*L3;
	H[3] = 0.5*(3 * L4 - 1)*(3 * L4 - 2)*L4;
	H[4] = 9.0 / 2.0*(3 * L1 - 1)*L1*L2;
	H[5] = 9.0 / 2.0*(3 * L2 - 1)*L1*L2;
	H[6] = 9.0 / 2.0*(3 * L2 - 1)*L2*L3;
	H[7] = 9.0 / 2.0*(3 * L3 - 1)*L2*L3;
	H[8] = 9.0 / 2.0*(3 * L1 - 1)*L1*L3;
	H[9] = 9.0 / 2.0*(3 * L3 - 1)*L1*L3;
	H[10] = 9.0 / 2.0*(3 * L1 - 1)*L1*L4;
	H[11] = 9.0 / 2.0*(3 * L4 - 1)*L1*L4;
	H[12] = 9.0 / 2.0*(3 * L2 - 1)*L2*L4;
	H[13] = 9.0 / 2.0*(3 * L4 - 1)*L2*L4;
	H[14] = 9.0 / 2.0*(3 * L3 - 1)*L3*L4;
	H[15] = 9.0 / 2.0*(3 * L4 - 1)*L3*L4;
	H[16] = 27.0*L1*L2*L4;
	H[17] = 27.0*L2*L3*L4;
	H[18] = 27.0*L1*L3*L4;
	H[19] = 27.0*L1*L2*L3;
}

//-----------------------------------------------------------------------------
//! values of shape function derivatives
void FETet20::shape_deriv(double* Hr, double* Hs, double* Ht, double r, double s, double t)
{
	double L1 = 1.0 - r - s - t;
	double L2 = r;
	double L3 = s;
	double L4 = t;

	Hr[0] = -3. / 2.*(3 * L1 - 2)*L1 - 3. / 2.*(3 * L1 - 1)*L1 - 0.5*(3 * L1 - 1)*(3 * L1 - 2);
	Hr[1] = 3. / 2.*(3 * L2 - 2)*L2 + 3. / 2.*(3 * L2 - 1)*L2 + 0.5*(3 * L2 - 1)*(3 * L2 - 2);
	Hr[2] = 0.0;
	Hr[3] = 0.0;
	Hr[4] = -27. / 2.*L1*L2 - 9.0 / 2.0*(3 * L1 - 1)*L2 + 9.0 / 2.0*(3 * L1 - 1)*L1;
	Hr[5] = 27. / 2.*L1*L2 - 9.0 / 2.0*(3 * L2 - 1)*L2 + 9.0 / 2.0*(3 * L2 - 1)*L1;
	Hr[6] = 27. / 2.*L2*L3 + 9.0 / 2.0*(3 * L2 - 1)*L3;
	Hr[7] = 9.0 / 2.0*(3 * L3 - 1)*L3;
	Hr[8] = -27. / 2.*L1*L3 - 9.0 / 2.0*(3 * L1 - 1)*L3;
	Hr[9] = -9.0 / 2.0*(3 * L3 - 1)*L3;
	Hr[10] = -27. / 2.*L1*L4 - 9.0 / 2.0*(3 * L1 - 1)*L4;
	Hr[11] = -9. / 2.*(3 * L4 - 1)*L4;
	Hr[12] = 27. / 2.*L2*L4 + 9. / 2.*(3 * L2 - 1)*L4;
	Hr[13] = 9. / 2.*(3 * L4 - 1)*L4;
	Hr[14] = 0.0;
	Hr[15] = 0.0;
	Hr[16] = -27 * L2*L4 + 27 * L1*L4;
	Hr[17] = 27 * L3*L4;
	Hr[18] = -27 * L3*L4;
	Hr[19] = -27 * L2*L3 + 27 * L1*L3;

	Hs[0] = -3. / 2.*(3 * L1 - 2)*L1 - 3. / 2.*(3 * L1 - 1)*L1 - 0.5*(3 * L1 - 1)*(3 * L1 - 2);
	Hs[1] = 0.0;
	Hs[2] = 3. / 2.*(3 * L3 - 2)*L3 + 3. / 2.*(3 * L3 - 1)*L3 + 0.5*(3 * L3 - 1)*(3 * L3 - 2);
	Hs[3] = 0.0;
	Hs[4] = -27. / 2.*L1*L2 - 9. / 2.*(3 * L1 - 1)*L2;
	Hs[5] = -9. / 2.*(3 * L2 - 1)*L2;
	Hs[6] = 9. / 2.*(3 * L2 - 1)*L2;
	Hs[7] = 27. / 2.*L2*L3 + 9. / 2.*(3 * L3 - 1)*L2;
	Hs[8] = -27. / 2.*L1*L3 - 9. / 2.*(3 * L1 - 1)*L3 + 9. / 2.*(3 * L1 - 1)*L1;
	Hs[9] = 27. / 2.*L1*L3 - 9. / 2.*(3 * L3 - 1)*L3 + 9. / 2.*(3 * L3 - 1)*L1;
	Hs[10] = -27. / 2.*L1*L4 - 9. / 2.*(3 * L1 - 1)*L4;
	Hs[11] = -9. / 2.*(3 * L4 - 1)*L4;
	Hs[12] = 0.0;
	Hs[13] = 0.0;
	Hs[14] = 27. / 2.*L3*L4 + 9. / 2.*(3 * L3 - 1)*L4;
	Hs[15] = 9. / 2.*(3 * L4 - 1)*L4;
	Hs[16] = -27 * L2*L4;
	Hs[17] = 27 * L2*L4;
	Hs[18] = -27 * L3*L4 + 27 * L1*L4;
	Hs[19] = -27 * L2*L3 + 27 * L1*L2;

	Ht[0] = -3. / 2.*(3 * L1 - 2)*L1 - 3. / 2.*(3 * L1 - 1)*L1 - 0.5*(3 * L1 - 1)*(3 * L1 - 2);
	Ht[1] = 0.0;
	Ht[2] = 0.0;
	Ht[3] = 3. / 2.*(3 * L4 - 2)*L4 + 3. / 2.*(3 * L4 - 1)*L4 + 0.5*(3 * L4 - 1)*(3 * L4 - 2);
	Ht[4] = -27. / 2.*L1*L2 - 9. / 2.*(3 * L1 - 1)*L2;
	Ht[5] = -9. / 2.*(3 * L2 - 1)*L2;
	Ht[6] = 0.0;
	Ht[7] = 0.0;
	Ht[8] = -27. / 2.*L1*L3 - 9. / 2.*(3 * L1 - 1)*L3;
	Ht[9] = -9. / 2.*(3 * L3 - 1)*L3;
	Ht[10] = -27. / 2.*L1*L4 - 9. / 2.*(3 * L1 - 1)*L4 + 9. / 2.*(3 * L1 - 1)*L1;
	Ht[11] = 27. / 2.*L1*L4 - 9. / 2.*(3 * L4 - 1)*L4 + 9. / 2.*(3 * L4 - 1)*L1;
	Ht[12] = 9. / 2.*(3 * L2 - 1)*L2;
	Ht[13] = 27. / 2.*L2*L4 + 9. / 2.*(3 * L4 - 1)*L2;
	Ht[14] = 9. / 2.*(3 * L3 - 1)*L3;
	Ht[15] = 27. / 2.*L3*L4 + 9. / 2.*(3 * L4 - 1)*L3;
	Ht[16] = -27 * L2*L4 + 27 * L1*L2;
	Ht[17] = 27 * L2*L3;
	Ht[18] = -27 * L3*L4 + 27 * L1*L3;
	Ht[19] = -27 * L2*L3;
}

//-----------------------------------------------------------------------------
//! \todo implement this
void FETet20::shape_deriv2(double* Hrr, double* Hss, double* Htt, double* Hrs, double* Hst, double* Hrt, double r, double s, double t)
{
	// do this
}


//=============================================================================
//              H E X 2 0
//=============================================================================

//-----------------------------------------------------------------------------
void FEHex20::shape_fnc(double* H, double r, double s, double t)
{
	H[8] = 0.25*(1 - r*r)*(1 - s)*(1 - t);
	H[9] = 0.25*(1 - s*s)*(1 + r)*(1 - t);
	H[10] = 0.25*(1 - r*r)*(1 + s)*(1 - t);
	H[11] = 0.25*(1 - s*s)*(1 - r)*(1 - t);
	H[12] = 0.25*(1 - r*r)*(1 - s)*(1 + t);
	H[13] = 0.25*(1 - s*s)*(1 + r)*(1 + t);
	H[14] = 0.25*(1 - r*r)*(1 + s)*(1 + t);
	H[15] = 0.25*(1 - s*s)*(1 - r)*(1 + t);
	H[16] = 0.25*(1 - t*t)*(1 - r)*(1 - s);
	H[17] = 0.25*(1 - t*t)*(1 + r)*(1 - s);
	H[18] = 0.25*(1 - t*t)*(1 + r)*(1 + s);
	H[19] = 0.25*(1 - t*t)*(1 - r)*(1 + s);

	H[0] = 0.125*(1 - r)*(1 - s)*(1 - t) - 0.5*(H[8] + H[11] + H[16]);
	H[1] = 0.125*(1 + r)*(1 - s)*(1 - t) - 0.5*(H[8] + H[9] + H[17]);
	H[2] = 0.125*(1 + r)*(1 + s)*(1 - t) - 0.5*(H[9] + H[10] + H[18]);
	H[3] = 0.125*(1 - r)*(1 + s)*(1 - t) - 0.5*(H[10] + H[11] + H[19]);
	H[4] = 0.125*(1 - r)*(1 - s)*(1 + t) - 0.5*(H[12] + H[15] + H[16]);
	H[5] = 0.125*(1 + r)*(1 - s)*(1 + t) - 0.5*(H[12] + H[13] + H[17]);
	H[6] = 0.125*(1 + r)*(1 + s)*(1 + t) - 0.5*(H[13] + H[14] + H[18]);
	H[7] = 0.125*(1 - r)*(1 + s)*(1 + t) - 0.5*(H[14] + H[15] + H[19]);
}

//-----------------------------------------------------------------------------
void FEHex20::shape_deriv(double* Hr, double* Hs, double* Ht, double r, double s, double t)
{
	Hr[8] = -0.5*r*(1 - s)*(1 - t);
	Hr[9] = 0.25*(1 - s*s)*(1 - t);
	Hr[10] = -0.5*r*(1 + s)*(1 - t);
	Hr[11] = -0.25*(1 - s*s)*(1 - t);
	Hr[12] = -0.5*r*(1 - s)*(1 + t);
	Hr[13] = 0.25*(1 - s*s)*(1 + t);
	Hr[14] = -0.5*r*(1 + s)*(1 + t);
	Hr[15] = -0.25*(1 - s*s)*(1 + t);
	Hr[16] = -0.25*(1 - t*t)*(1 - s);
	Hr[17] = 0.25*(1 - t*t)*(1 - s);
	Hr[18] = 0.25*(1 - t*t)*(1 + s);
	Hr[19] = -0.25*(1 - t*t)*(1 + s);

	Hr[0] = -0.125*(1 - s)*(1 - t) - 0.5*(Hr[8] + Hr[11] + Hr[16]);
	Hr[1] = 0.125*(1 - s)*(1 - t) - 0.5*(Hr[8] + Hr[9] + Hr[17]);
	Hr[2] = 0.125*(1 + s)*(1 - t) - 0.5*(Hr[9] + Hr[10] + Hr[18]);
	Hr[3] = -0.125*(1 + s)*(1 - t) - 0.5*(Hr[10] + Hr[11] + Hr[19]);
	Hr[4] = -0.125*(1 - s)*(1 + t) - 0.5*(Hr[12] + Hr[15] + Hr[16]);
	Hr[5] = 0.125*(1 - s)*(1 + t) - 0.5*(Hr[12] + Hr[13] + Hr[17]);
	Hr[6] = 0.125*(1 + s)*(1 + t) - 0.5*(Hr[13] + Hr[14] + Hr[18]);
	Hr[7] = -0.125*(1 + s)*(1 + t) - 0.5*(Hr[14] + Hr[15] + Hr[19]);

	Hs[8] = -0.25*(1 - r*r)*(1 - t);
	Hs[9] = -0.5*s*(1 + r)*(1 - t);
	Hs[10] = 0.25*(1 - r*r)*(1 - t);
	Hs[11] = -0.5*s*(1 - r)*(1 - t);
	Hs[12] = -0.25*(1 - r*r)*(1 + t);
	Hs[13] = -0.5*s*(1 + r)*(1 + t);
	Hs[14] = 0.25*(1 - r*r)*(1 + t);
	Hs[15] = -0.5*s*(1 - r)*(1 + t);
	Hs[16] = -0.25*(1 - t*t)*(1 - r);
	Hs[17] = -0.25*(1 - t*t)*(1 + r);
	Hs[18] = 0.25*(1 - t*t)*(1 + r);
	Hs[19] = 0.25*(1 - t*t)*(1 - r);

	Hs[0] = -0.125*(1 - r)*(1 - t) - 0.5*(Hs[8] + Hs[11] + Hs[16]);
	Hs[1] = -0.125*(1 + r)*(1 - t) - 0.5*(Hs[8] + Hs[9] + Hs[17]);
	Hs[2] = 0.125*(1 + r)*(1 - t) - 0.5*(Hs[9] + Hs[10] + Hs[18]);
	Hs[3] = 0.125*(1 - r)*(1 - t) - 0.5*(Hs[10] + Hs[11] + Hs[19]);
	Hs[4] = -0.125*(1 - r)*(1 + t) - 0.5*(Hs[12] + Hs[15] + Hs[16]);
	Hs[5] = -0.125*(1 + r)*(1 + t) - 0.5*(Hs[12] + Hs[13] + Hs[17]);
	Hs[6] = 0.125*(1 + r)*(1 + t) - 0.5*(Hs[13] + Hs[14] + Hs[18]);
	Hs[7] = 0.125*(1 - r)*(1 + t) - 0.5*(Hs[14] + Hs[15] + Hs[19]);

	Ht[8] = -0.25*(1 - r*r)*(1 - s);
	Ht[9] = -0.25*(1 - s*s)*(1 + r);
	Ht[10] = -0.25*(1 - r*r)*(1 + s);
	Ht[11] = -0.25*(1 - s*s)*(1 - r);
	Ht[12] = 0.25*(1 - r*r)*(1 - s);
	Ht[13] = 0.25*(1 - s*s)*(1 + r);
	Ht[14] = 0.25*(1 - r*r)*(1 + s);
	Ht[15] = 0.25*(1 - s*s)*(1 - r);
	Ht[16] = -0.5*t*(1 - r)*(1 - s);
	Ht[17] = -0.5*t*(1 + r)*(1 - s);
	Ht[18] = -0.5*t*(1 + r)*(1 + s);
	Ht[19] = -0.5*t*(1 - r)*(1 + s);

	Ht[0] = -0.125*(1 - r)*(1 - s) - 0.5*(Ht[8] + Ht[11] + Ht[16]);
	Ht[1] = -0.125*(1 + r)*(1 - s) - 0.5*(Ht[8] + Ht[9] + Ht[17]);
	Ht[2] = -0.125*(1 + r)*(1 + s) - 0.5*(Ht[9] + Ht[10] + Ht[18]);
	Ht[3] = -0.125*(1 - r)*(1 + s) - 0.5*(Ht[10] + Ht[11] + Ht[19]);
	Ht[4] = 0.125*(1 - r)*(1 - s) - 0.5*(Ht[12] + Ht[15] + Ht[16]);
	Ht[5] = 0.125*(1 + r)*(1 - s) - 0.5*(Ht[12] + Ht[13] + Ht[17]);
	Ht[6] = 0.125*(1 + r)*(1 + s) - 0.5*(Ht[13] + Ht[14] + Ht[18]);
	Ht[7] = 0.125*(1 - r)*(1 + s) - 0.5*(Ht[14] + Ht[15] + Ht[19]);
}

//-----------------------------------------------------------------------------
void FEHex20::shape_deriv2(double* Hrr, double* Hss, double* Htt, double* Hrs, double* Hst, double* Hrt, double r, double s, double t)
{
	// Hrr
	Hrr[8] = -0.5*(1 - s)*(1 - t);
	Hrr[9] = 0.;
	Hrr[10] = -0.5*(1 + s)*(1 - t);
	Hrr[11] = 0.;
	Hrr[12] = -0.5*(1 - s)*(1 + t);
	Hrr[13] = 0.;
	Hrr[14] = -0.5*(1 + s)*(1 + t);
	Hrr[15] = 0.;
	Hrr[16] = 0.;
	Hrr[17] = 0.;
	Hrr[18] = 0.;
	Hrr[19] = 0.;

	Hrr[0] = -0.5*(Hrr[8] + Hrr[11] + Hrr[16]);
	Hrr[1] = -0.5*(Hrr[8] + Hrr[9] + Hrr[17]);
	Hrr[2] = -0.5*(Hrr[9] + Hrr[10] + Hrr[18]);
	Hrr[3] = -0.5*(Hrr[10] + Hrr[11] + Hrr[19]);
	Hrr[4] = -0.5*(Hrr[12] + Hrr[15] + Hrr[16]);
	Hrr[5] = -0.5*(Hrr[12] + Hrr[13] + Hrr[17]);
	Hrr[6] = -0.5*(Hrr[13] + Hrr[14] + Hrr[18]);
	Hrr[7] = -0.5*(Hrr[14] + Hrr[15] + Hrr[19]);

	// Hss
	Hss[8] = 0.;
	Hss[9] = -0.5*(1 + r)*(1 - t);
	Hss[10] = 0.;
	Hss[11] = -0.5*(1 - r)*(1 - t);
	Hss[12] = 0.;
	Hss[13] = -0.5*(1 + r)*(1 + t);
	Hss[14] = 0.;
	Hss[15] = -0.5*(1 - r)*(1 + t);
	Hss[16] = 0.;
	Hss[17] = 0.;
	Hss[18] = 0.;
	Hss[19] = 0.;

	Hss[0] = -0.5*(Hss[8] + Hss[11] + Hss[16]);
	Hss[1] = -0.5*(Hss[8] + Hss[9] + Hss[17]);
	Hss[2] = -0.5*(Hss[9] + Hss[10] + Hss[18]);
	Hss[3] = -0.5*(Hss[10] + Hss[11] + Hss[19]);
	Hss[4] = -0.5*(Hss[12] + Hss[15] + Hss[16]);
	Hss[5] = -0.5*(Hss[12] + Hss[13] + Hss[17]);
	Hss[6] = -0.5*(Hss[13] + Hss[14] + Hss[18]);
	Hss[7] = -0.5*(Hss[14] + Hss[15] + Hss[19]);

	// Htt
	Htt[8] = 0.;
	Htt[9] = 0.;
	Htt[10] = 0.;
	Htt[11] = 0.;
	Htt[12] = 0.;
	Htt[13] = 0.;
	Htt[14] = 0.;
	Htt[15] = 0.;
	Htt[16] = -0.5*(1 - r)*(1 - s);
	Htt[17] = -0.5*(1 + r)*(1 - s);
	Htt[18] = -0.5*(1 + r)*(1 + s);
	Htt[19] = -0.5*(1 - r)*(1 + s);

	Htt[0] = -0.5*(Htt[8] + Htt[11] + Htt[16]);
	Htt[1] = -0.5*(Htt[8] + Htt[9] + Htt[17]);
	Htt[2] = -0.5*(Htt[9] + Htt[10] + Htt[18]);
	Htt[3] = -0.5*(Htt[10] + Htt[11] + Htt[19]);
	Htt[4] = -0.5*(Htt[12] + Htt[15] + Htt[16]);
	Htt[5] = -0.5*(Htt[12] + Htt[13] + Htt[17]);
	Htt[6] = -0.5*(Htt[13] + Htt[14] + Htt[18]);
	Htt[7] = -0.5*(Htt[14] + Htt[15] + Htt[19]);

	// Hrs
	Hrs[8] = 0.5*r*(1 - t);
	Hrs[9] = -0.5*s*(1 - t);
	Hrs[10] = -0.5*r*(1 - t);
	Hrs[11] = 0.5*s*(1 - t);
	Hrs[12] = 0.5*r*(1 + t);
	Hrs[13] = -0.5*s*(1 + t);
	Hrs[14] = -0.5*r*(1 + t);
	Hrs[15] = 0.5*s*(1 + t);
	Hrs[16] = 0.25*(1 - t*t);
	Hrs[17] = -0.25*(1 - t*t);
	Hrs[18] = 0.25*(1 - t*t);
	Hrs[19] = -0.25*(1 - t*t);

	Hrs[0] = 0.125*(1 - t) - 0.5*(Hrs[8] + Hrs[11] + Hrs[16]);
	Hrs[1] = -0.125*(1 - t) - 0.5*(Hrs[8] + Hrs[9] + Hrs[17]);
	Hrs[2] = 0.125*(1 - t) - 0.5*(Hrs[9] + Hrs[10] + Hrs[18]);
	Hrs[3] = -0.125*(1 - t) - 0.5*(Hrs[10] + Hrs[11] + Hrs[19]);
	Hrs[4] = 0.125*(1 + t) - 0.5*(Hrs[12] + Hrs[15] + Hrs[16]);
	Hrs[5] = -0.125*(1 + t) - 0.5*(Hrs[12] + Hrs[13] + Hrs[17]);
	Hrs[6] = 0.125*(1 + t) - 0.5*(Hrs[13] + Hrs[14] + Hrs[18]);
	Hrs[7] = -0.125*(1 + t) - 0.5*(Hrs[14] + Hrs[15] + Hrs[19]);

	// Hst
	Hst[8] = 0.25*(1 - r*r);
	Hst[9] = 0.5*s*(1 + r);
	Hst[10] = -0.25*(1 - r*r);
	Hst[11] = 0.5*s*(1 - r);
	Hst[12] = -0.25*(1 - r*r);
	Hst[13] = -0.5*s*(1 + r);
	Hst[14] = 0.25*(1 - r*r);
	Hst[15] = -0.5*s*(1 - r);
	Hst[16] = 0.5*t*(1 - r);
	Hst[17] = 0.5*t*(1 + r);
	Hst[18] = -0.5*t*(1 + r);
	Hst[19] = -0.5*t*(1 - r);

	Hst[0] = 0.125*(1 - r) - 0.5*(Hst[8] + Hst[11] + Hst[16]);
	Hst[1] = 0.125*(1 + r) - 0.5*(Hst[8] + Hst[9] + Hst[17]);
	Hst[2] = -0.125*(1 + r) - 0.5*(Hst[9] + Hst[10] + Hst[18]);
	Hst[3] = -0.125*(1 - r) - 0.5*(Hst[10] + Hst[11] + Hst[19]);
	Hst[4] = -0.125*(1 - r) - 0.5*(Hst[12] + Hst[15] + Hst[16]);
	Hst[5] = -0.125*(1 + r) - 0.5*(Hst[12] + Hst[13] + Hst[17]);
	Hst[6] = 0.125*(1 + r) - 0.5*(Hst[13] + Hst[14] + Hst[18]);
	Hst[7] = 0.125*(1 - r) - 0.5*(Hst[14] + Hst[15] + Hst[19]);

	// Hrt
	Hrt[8] = 0.5*r*(1 - s);
	Hrt[9] = -0.25*(1 - s*s);
	Hrt[10] = 0.5*r*(1 + s);
	Hrt[11] = 0.25*(1 - s*s);
	Hrt[12] = -0.5*r*(1 - s);
	Hrt[13] = 0.25*(1 - s*s);
	Hrt[14] = -0.5*r*(1 + s);
	Hrt[15] = -0.25*(1 - s*s);
	Hrt[16] = 0.5*t*(1 - s);
	Hrt[17] = -0.5*t*(1 - s);
	Hrt[18] = -0.5*t*(1 + s);
	Hrt[19] = 0.5*t*(1 + s);

	Hrt[0] = 0.125*(1 - s) - 0.5*(Hrt[8] + Hrt[11] + Hrt[16]);
	Hrt[1] = -0.125*(1 - s) - 0.5*(Hrt[8] + Hrt[9] + Hrt[17]);
	Hrt[2] = -0.125*(1 + s) - 0.5*(Hrt[9] + Hrt[10] + Hrt[18]);
	Hrt[3] = 0.125*(1 + s) - 0.5*(Hrt[10] + Hrt[11] + Hrt[19]);
	Hrt[4] = -0.125*(1 - s) - 0.5*(Hrt[12] + Hrt[15] + Hrt[16]);
	Hrt[5] = 0.125*(1 - s) - 0.5*(Hrt[12] + Hrt[13] + Hrt[17]);
	Hrt[6] = 0.125*(1 + s) - 0.5*(Hrt[13] + Hrt[14] + Hrt[18]);
	Hrt[7] = -0.125*(1 + s) - 0.5*(Hrt[14] + Hrt[15] + Hrt[19]);
}

//=============================================================================
//              H E X 2 7
//=============================================================================
// Lookup table for 27-node hex that maps a node index into a triplet that identifies
// the shape functions.
static int HEX27_LUT[27][3] = {
	{ 0,0,0 },
	{ 1,0,0 },
	{ 1,1,0 },
	{ 0,1,0 },
	{ 0,0,1 },
	{ 1,0,1 },
	{ 1,1,1 },
	{ 0,1,1 },
	{ 2,0,0 },
	{ 1,2,0 },
	{ 2,1,0 },
	{ 0,2,0 },
	{ 2,0,1 },
	{ 1,2,1 },
	{ 2,1,1 },
	{ 0,2,1 },
	{ 0,0,2 },
	{ 1,0,2 },
	{ 1,1,2 },
	{ 0,1,2 },
	{ 2,0,2 },
	{ 1,2,2 },
	{ 2,1,2 },
	{ 0,2,2 },
	{ 2,2,0 },
	{ 2,2,1 },
	{ 2,2,2 }
};

//-----------------------------------------------------------------------------
//! values of shape functions
void FEHex27::shape_fnc(double* H, double r, double s, double t)
{
	double R[3] = { 0.5*r*(r - 1.0), 0.5*r*(r + 1.0), 1.0 - r*r };
	double S[3] = { 0.5*s*(s - 1.0), 0.5*s*(s + 1.0), 1.0 - s*s };
	double T[3] = { 0.5*t*(t - 1.0), 0.5*t*(t + 1.0), 1.0 - t*t };

	H[0] = R[0] * S[0] * T[0];
	H[1] = R[1] * S[0] * T[0];
	H[2] = R[1] * S[1] * T[0];
	H[3] = R[0] * S[1] * T[0];
	H[4] = R[0] * S[0] * T[1];
	H[5] = R[1] * S[0] * T[1];
	H[6] = R[1] * S[1] * T[1];
	H[7] = R[0] * S[1] * T[1];
	H[8] = R[2] * S[0] * T[0];
	H[9] = R[1] * S[2] * T[0];
	H[10] = R[2] * S[1] * T[0];
	H[11] = R[0] * S[2] * T[0];
	H[12] = R[2] * S[0] * T[1];
	H[13] = R[1] * S[2] * T[1];
	H[14] = R[2] * S[1] * T[1];
	H[15] = R[0] * S[2] * T[1];
	H[16] = R[0] * S[0] * T[2];
	H[17] = R[1] * S[0] * T[2];
	H[18] = R[1] * S[1] * T[2];
	H[19] = R[0] * S[1] * T[2];
	H[20] = R[2] * S[0] * T[2];
	H[21] = R[1] * S[2] * T[2];
	H[22] = R[2] * S[1] * T[2];
	H[23] = R[0] * S[2] * T[2];
	H[24] = R[2] * S[2] * T[0];
	H[25] = R[2] * S[2] * T[1];
	H[26] = R[2] * S[2] * T[2];
}

//-----------------------------------------------------------------------------
//! values of shape function derivatives
void FEHex27::shape_deriv(double* Hr, double* Hs, double* Ht, double r, double s, double t)
{
	double R[3] = { 0.5*r*(r - 1.0), 0.5*r*(r + 1.0), 1.0 - r*r };
	double S[3] = { 0.5*s*(s - 1.0), 0.5*s*(s + 1.0), 1.0 - s*s };
	double T[3] = { 0.5*t*(t - 1.0), 0.5*t*(t + 1.0), 1.0 - t*t };

	double DR[3] = { r - 0.5, r + 0.5, -2.0*r };
	double DS[3] = { s - 0.5, s + 0.5, -2.0*s };
	double DT[3] = { t - 0.5, t + 0.5, -2.0*t };

	Hr[0] = DR[0] * S[0] * T[0]; Hs[0] = R[0] * DS[0] * T[0]; Ht[0] = R[0] * S[0] * DT[0];
	Hr[1] = DR[1] * S[0] * T[0]; Hs[1] = R[1] * DS[0] * T[0];	Ht[1] = R[1] * S[0] * DT[0];
	Hr[2] = DR[1] * S[1] * T[0]; Hs[2] = R[1] * DS[1] * T[0];	Ht[2] = R[1] * S[1] * DT[0];
	Hr[3] = DR[0] * S[1] * T[0]; Hs[3] = R[0] * DS[1] * T[0];	Ht[3] = R[0] * S[1] * DT[0];
	Hr[4] = DR[0] * S[0] * T[1]; Hs[4] = R[0] * DS[0] * T[1];	Ht[4] = R[0] * S[0] * DT[1];
	Hr[5] = DR[1] * S[0] * T[1]; Hs[5] = R[1] * DS[0] * T[1];	Ht[5] = R[1] * S[0] * DT[1];
	Hr[6] = DR[1] * S[1] * T[1]; Hs[6] = R[1] * DS[1] * T[1];	Ht[6] = R[1] * S[1] * DT[1];
	Hr[7] = DR[0] * S[1] * T[1]; Hs[7] = R[0] * DS[1] * T[1];	Ht[7] = R[0] * S[1] * DT[1];
	Hr[8] = DR[2] * S[0] * T[0]; Hs[8] = R[2] * DS[0] * T[0];	Ht[8] = R[2] * S[0] * DT[0];
	Hr[9] = DR[1] * S[2] * T[0]; Hs[9] = R[1] * DS[2] * T[0];	Ht[9] = R[1] * S[2] * DT[0];
	Hr[10] = DR[2] * S[1] * T[0]; Hs[10] = R[2] * DS[1] * T[0];	Ht[10] = R[2] * S[1] * DT[0];
	Hr[11] = DR[0] * S[2] * T[0]; Hs[11] = R[0] * DS[2] * T[0];	Ht[11] = R[0] * S[2] * DT[0];
	Hr[12] = DR[2] * S[0] * T[1]; Hs[12] = R[2] * DS[0] * T[1];	Ht[12] = R[2] * S[0] * DT[1];
	Hr[13] = DR[1] * S[2] * T[1]; Hs[13] = R[1] * DS[2] * T[1];	Ht[13] = R[1] * S[2] * DT[1];
	Hr[14] = DR[2] * S[1] * T[1]; Hs[14] = R[2] * DS[1] * T[1];	Ht[14] = R[2] * S[1] * DT[1];
	Hr[15] = DR[0] * S[2] * T[1]; Hs[15] = R[0] * DS[2] * T[1];	Ht[15] = R[0] * S[2] * DT[1];
	Hr[16] = DR[0] * S[0] * T[2]; Hs[16] = R[0] * DS[0] * T[2];	Ht[16] = R[0] * S[0] * DT[2];
	Hr[17] = DR[1] * S[0] * T[2]; Hs[17] = R[1] * DS[0] * T[2];	Ht[17] = R[1] * S[0] * DT[2];
	Hr[18] = DR[1] * S[1] * T[2]; Hs[18] = R[1] * DS[1] * T[2];	Ht[18] = R[1] * S[1] * DT[2];
	Hr[19] = DR[0] * S[1] * T[2]; Hs[19] = R[0] * DS[1] * T[2];	Ht[19] = R[0] * S[1] * DT[2];
	Hr[20] = DR[2] * S[0] * T[2]; Hs[20] = R[2] * DS[0] * T[2];	Ht[20] = R[2] * S[0] * DT[2];
	Hr[21] = DR[1] * S[2] * T[2]; Hs[21] = R[1] * DS[2] * T[2];	Ht[21] = R[1] * S[2] * DT[2];
	Hr[22] = DR[2] * S[1] * T[2]; Hs[22] = R[2] * DS[1] * T[2];	Ht[22] = R[2] * S[1] * DT[2];
	Hr[23] = DR[0] * S[2] * T[2]; Hs[23] = R[0] * DS[2] * T[2];	Ht[23] = R[0] * S[2] * DT[2];
	Hr[24] = DR[2] * S[2] * T[0]; Hs[24] = R[2] * DS[2] * T[0];	Ht[24] = R[2] * S[2] * DT[0];
	Hr[25] = DR[2] * S[2] * T[1]; Hs[25] = R[2] * DS[2] * T[1];	Ht[25] = R[2] * S[2] * DT[1];
	Hr[26] = DR[2] * S[2] * T[2]; Hs[26] = R[2] * DS[2] * T[2];	Ht[26] = R[2] * S[2] * DT[2];
}

//-----------------------------------------------------------------------------
//! values of shape function second derivatives
void FEHex27::shape_deriv2(double* Hrr, double* Hss, double* Htt, double* Hrs, double* Hst, double* Hrt, double r, double s, double t)
{
	double NR[3] = { 0.5*r*(r - 1.0), 0.5*r*(r + 1.0), 1.0 - r*r };
	double NS[3] = { 0.5*s*(s - 1.0), 0.5*s*(s + 1.0), 1.0 - s*s };
	double NT[3] = { 0.5*t*(t - 1.0), 0.5*t*(t + 1.0), 1.0 - t*t };

	double DR[3] = { r - 0.5, r + 0.5, -2.0*r };
	double DS[3] = { s - 0.5, s + 0.5, -2.0*s };
	double DT[3] = { t - 0.5, t + 0.5, -2.0*t };

	double HR[3] = { 1.0, 1.0, -2.0 };
	double HS[3] = { 1.0, 1.0, -2.0 };
	double HT[3] = { 1.0, 1.0, -2.0 };

	for (int a = 0; a<27; ++a)
	{
		int i = HEX27_LUT[a][0];
		int j = HEX27_LUT[a][1];
		int k = HEX27_LUT[a][2];

		Hrr[a] = HR[i] * NS[j] * NT[k];
		Hss[a] = NR[i] * HS[j] * NT[k];
		Htt[a] = NR[i] * NS[j] * HT[k];

		Hrs[a] = DR[i] * DS[j] * NT[k];
		Hst[a] = NR[i] * DS[j] * DT[k];
		Hrt[a] = DR[i] * NS[j] * DT[k];
	}
}

//=============================================================================
//              P Y R A 5
//=============================================================================

//! values of shape functions
void FEPyra5::shape_fnc(double* H, double r, double s, double t)
{
	H[0] = 0.125*(1.0 - r)*(1.0 - s)*(1.0 - t);
	H[1] = 0.125*(1.0 + r)*(1.0 - s)*(1.0 - t);
	H[2] = 0.125*(1.0 + r)*(1.0 + s)*(1.0 - t);
	H[3] = 0.125*(1.0 - r)*(1.0 + s)*(1.0 - t);
	H[4] = 0.5*(1.0 + t);
}

//! values of shape function derivatives
void FEPyra5::shape_deriv(double* Hr, double* Hs, double* Ht, double r, double s, double t)
{
	Hr[0] = -0.125*(1.0 - s)*(1.0 - t);
	Hr[1] = 0.125*(1.0 - s)*(1.0 - t);
	Hr[2] = 0.125*(1.0 + s)*(1.0 - t);
	Hr[3] = -0.125*(1.0 + s)*(1.0 - t);
	Hr[4] = 0.0;

	Hs[0] = -0.125*(1.0 - r)*(1.0 - t);
	Hs[1] = -0.125*(1.0 + r)*(1.0 - t);
	Hs[2] = 0.125*(1.0 + r)*(1.0 - t);
	Hs[3] = 0.125*(1.0 - r)*(1.0 - t);
	Hs[4] = 0.0;

	Ht[0] = -0.125*(1.0 - r)*(1.0 - s);
	Ht[1] = -0.125*(1.0 + r)*(1.0 - s);
	Ht[2] = -0.125*(1.0 + r)*(1.0 + s);
	Ht[3] = -0.125*(1.0 - r)*(1.0 + s);
	Ht[4] = 0.5;
}

//! values of shape function second derivatives
void FEPyra5::shape_deriv2(double* Hrr, double* Hss, double* Htt, double* Hrs, double* Hst, double* Hrt, double r, double s, double t)
{
	Hrr[0] = 0.0; Hss[0] = 0.0; Htt[0] = 0.0;
	Hrr[1] = 0.0; Hss[1] = 0.0; Htt[1] = 0.0;
	Hrr[2] = 0.0; Hss[2] = 0.0; Htt[2] = 0.0;
	Hrr[3] = 0.0; Hss[3] = 0.0; Htt[3] = 0.0;
	Hrr[4] = 0.0; Hss[4] = 0.0; Htt[4] = 0.0;

	Hrs[0] = 0.125*(1.0 - t); Hrt[0] = 0.125*(1.0 - s); Hst[0] = 0.125*(1.0 - r);
	Hrs[1] = -0.125*(1.0 - t); Hrt[1] = -0.125*(1.0 - s); Hst[1] = 0.125*(1.0 + r);
	Hrs[2] = 0.125*(1.0 - t); Hrt[2] = -0.125*(1.0 + s); Hst[2] = -0.125*(1.0 + r);
	Hrs[3] = -0.125*(1.0 - t); Hrt[3] = 0.125*(1.0 + s); Hst[3] = -0.125*(1.0 - r);
	Hrs[4] = 0.0; Hrt[4] = 0.0; Hst[4] = 0.0;
}

//=============================================================================
//              P Y R A 1 3
//=============================================================================

//! values of shape functions
void FEPyra13::shape_fnc(double* H, double r, double s, double t)
{
    H[5] = 0.25*(1 - r*r)*(1 - s)*(1 - t);
    H[6] = 0.25*(1 - s*s)*(1 + r)*(1 - t);
    H[7] = 0.25*(1 - r*r)*(1 + s)*(1 - t);
    H[8] = 0.25*(1 - s*s)*(1 - r)*(1 - t);
    H[9] = 0.25*(1 - t*t)*(1 - r)*(1 - s);
    H[10] = 0.25*(1 - t*t)*(1 + r)*(1 - s);
    H[11] = 0.25*(1 - t*t)*(1 + r)*(1 + s);
    H[12] = 0.25*(1 - t*t)*(1 - r)*(1 + s);
    
    H[0] = 0.125*(1 - r)*(1 - s)*(1 - t) - 0.5*(H[5] + H[8] + H[9]);
    H[1] = 0.125*(1 + r)*(1 - s)*(1 - t) - 0.5*(H[5] + H[6] + H[10]);
    H[2] = 0.125*(1 + r)*(1 + s)*(1 - t) - 0.5*(H[6] + H[7] + H[11]);
    H[3] = 0.125*(1 - r)*(1 + s)*(1 - t) - 0.5*(H[7] + H[8] + H[12]);
    H[4] = 0.5*t*(1 + t);
}

//! values of shape function derivatives
void FEPyra13::shape_deriv(double* Hr, double* Hs, double* Ht, double r, double s, double t)
{
    Hr[ 0] = 0.125 + r*(0.25 + s*(-0.25 + 0.25*t) - 0.25*t) + s*s*(-0.125 + 0.125*t) +
    s*(-0.125 + 0.125*t)*t - 0.125*t*t;
    Hr[ 1] = -0.125 + r*(0.25 + s*(-0.25 + 0.25*t) - 0.25*t) + s*s*(0.125 - 0.125*t) +
    s*(0.125 - 0.125*t)*t + 0.125*t*t;
    Hr[ 2] = -0.125 + r*(0.25 + s*(0.25 - 0.25*t) - 0.25*t) + s*s*(0.125 - 0.125*t) +
    s*(-0.125 + 0.125*t)*t + 0.125*t*t;
    Hr[ 3] = 0.125 + r*(0.25 + s*(0.25 - 0.25*t) - 0.25*t) + s*s*(-0.125 + 0.125*t) +
    s*(0.125 - 0.125*t)*t - 0.125*t*t;
    Hr[ 4] = 0;
    Hr[ 5] = -0.5*r*(-1. + s)*(-1. + t);
    Hr[ 6] = 0.25*(-1. + s*s)*(-1. + t);
    Hr[ 7] = 0.5*r*(1. + s)*(-1. + t);
    Hr[ 8] = -0.25*(-1. + s*s)*(-1. + t);
    Hr[ 9] = -0.25*(-1. + s)*(-1. + t*t);
    Hr[10] = 0.25*(-1. + s)*(-1. + t*t);
    Hr[11] = -0.25*(1. + s)*(-1. + t*t);
    Hr[12] = 0.25*(1. + s)*(-1. + t*t);
    
    Hs[ 0] = 0.125 + s*(0.25 - 0.25*t) + r*r*(-0.125 + 0.125*t) - 0.125*t*t +
    r*(s*(-0.25 + 0.25*t) + (-0.125 + 0.125*t)*t);
    Hs[ 1] = 0.125 + s*(0.25 - 0.25*t) + r*r*(-0.125 + 0.125*t) - 0.125*t*t +
    r*(s*(0.25 - 0.25*t) + (0.125 - 0.125*t)*t);
    Hs[ 2] = -0.125 + s*(0.25 - 0.25*t) + r*r*(0.125 - 0.125*t) + 0.125*t*t +
    r*(s*(0.25 - 0.25*t) + (-0.125 + 0.125*t)*t);
    Hs[ 3] = -0.125 + s*(0.25 - 0.25*t) + r*r*(0.125 - 0.125*t) + 0.125*t*t +
    r*(s*(-0.25 + 0.25*t) + (0.125 - 0.125*t)*t);
    Hs[ 4] = 0;
    Hs[ 5] = -0.25*(-1. + r*r)*(-1. + t);
    Hs[ 6] = 0.5*(1. + r)*s*(-1. + t);
    Hs[ 7] = 0.25*(-1. + r*r)*(-1. + t);
    Hs[ 8] = -0.5*(-1. + r)*s*(-1. + t);
    Hs[ 9] = -0.25*(-1. + r)*(-1. + t*t);
    Hs[10] = 0.25*(1. + r)*(-1. + t*t);
    Hs[11] = -0.25*(1. + r)*(-1. + t*t);
    Hs[12] = 0.25*(-1. + r)*(-1. + t*t);

    Ht[ 0] = -0.125*(-1. + r)*(-1. + s) + 0.125*(-1. + r*r)*(-1. + s) +
    0.125*(-1. + r)*(-1. + s*s) + 0.25*(-1. + r)*(-1. + s)*t;
    Ht[ 1] = 0.125*(1. + r)*(-1. + s) + 0.125*(-1. + r*r)*(-1. + s) -
    0.125*(1. + r)*(-1. + s*s) - 0.25*(1. + r)*(-1. + s)*t;
    Ht[ 2] = -0.125*(1. + r)*(1. + s) - 0.125*(-1. + r*r)*(1. + s) -
    0.125*(1. + r)*(-1. + s*s) + 0.25*(1. + r)*(1. + s)*t;
    Ht[ 3] = 0.125*(-1. + r)*(1. + s) - 0.125*(-1. + r*r)*(1. + s) +
    0.125*(-1. + r)*(-1. + s*s) - 0.25*(-1. + r)*(1. + s)*t;
    Ht[ 4] = 0.5 + 1.*t;
    Ht[ 5] = -0.25*(-1. + r*r)*(-1. + s);
    Ht[ 6] = 0.25*(1. + r)*(-1. + s*s);
    Ht[ 7] = 0.25*(-1. + r*r)*(1. + s);
    Ht[ 8] = -0.25*(-1. + r)*(-1. + s*s);
    Ht[ 9] = -0.5*(-1. + r)*(-1. + s)*t;
    Ht[10] = 0.5*(1. + r)*(-1. + s)*t;
    Ht[11] = -0.5*(1. + r)*(1. + s)*t;
    Ht[12] = 0.5*(-1. + r)*(1. + s)*t;
}

//! values of shape function second derivatives
void FEPyra13::shape_deriv2(double* Hrr, double* Hss, double* Htt, double* Hrs, double* Hst, double* Hrt, double r, double s, double t)
{
    Hrr[ 0] = 0.25*(-1. + s)*(-1. + t);
    Hrr[ 1] = 0.25*(-1. + s)*(-1. + t);
    Hrr[ 2] = -0.25*(1. + s)*(-1. + t);
    Hrr[ 3] = -0.25*(1. + s)*(-1. + t);
    Hrr[ 4] = 0;
    Hrr[ 5] = -0.5*(-1. + s)*(-1. + t);
    Hrr[ 6] = 0.;
    Hrr[ 7] = 0.5*(1. + s)*(-1. + t);
    Hrr[ 8] = 0.;
    Hrr[ 9] = 0.;
    Hrr[10] = 0.;
    Hrr[11] = 0.;
    Hrr[12] = 0.;

    Hss[ 0] = 0.25*(-1. + r)*(-1. + t);
    Hss[ 1] = -0.25*(1. + r)*(-1. + t);
    Hss[ 2] = -0.25*(1. + r)*(-1. + t);
    Hss[ 3] = 0.25*(-1. + r)*(-1. + t);
    Hss[ 4] = 0;
    Hss[ 5] = 0;
    Hss[ 6] = 0.5*(1. + r)*(-1. + t);
    Hss[ 7] = 0;
    Hss[ 8] = -0.5*(-1. + r)*(-1. + t);
    Hss[ 9] = 0;
    Hss[10] = 0;
    Hss[11] = 0;
    Hss[12] = 0;

    Htt[ 0] = 0.25*(-1. + r)*(-1. + s);
    Htt[ 1] = -0.25*(1. + r)*(-1. + s);
    Htt[ 2] = 0.25*(1. + r)*(1. + s);
    Htt[ 3] = -0.25*(-1. + r)*(1. + s);
    Htt[ 4] = 1.;
    Htt[ 5] = 0;
    Htt[ 6] = 0;
    Htt[ 7] = 0;
    Htt[ 8] = 0;
    Htt[ 9] = -0.5*(-1. + r)*(-1. + s);
    Htt[10] = 0.5*(1. + r)*(-1. + s);
    Htt[11] = -0.5*(1. + r)*(1. + s);
    Htt[12] = 0.5*(-1. + r)*(1. + s);

    Hrs[ 0] = r*(-0.25 + 0.25*t) + s*(-0.25 + 0.25*t) - 0.125*t + 0.125*t*t;
    Hrs[ 1] = s*(0.25 - 0.25*t) + r*(-0.25 + 0.25*t) + 0.125*t - 0.125*t*t;
    Hrs[ 2] = r*(0.25 - 0.25*t) + s*(0.25 - 0.25*t) - 0.125*t + 0.125*t*t;
    Hrs[ 3] = r*(0.25 - 0.25*t) + s*(-0.25 + 0.25*t) + 0.125*t - 0.125*t*t;
    Hrs[ 4] = 0;
    Hrs[ 5] = -0.5*r*(-1. + t);
    Hrs[ 6] = 0.5*s*(-1. + t);
    Hrs[ 7] = 0.5*r*(-1. + t);
    Hrs[ 8] = -0.5*s*(-1. + t);
    Hrs[ 9] = -0.25*(-1. + t*t);
    Hrs[10] = 0.25*(-1. + t*t);
    Hrs[11] = -0.25*(-1. + t*t);
    Hrs[12] = 0.25*(-1. + t*t);
    
    Hrs[ 0] = 0.125*r*r - 0.25*s + r*(-0.125 + 0.25*s + 0.25*t) - 0.25*t;
    Hst[ 1] = 0.125*r*r - 0.25*s + r*(0.125 - 0.25*s - 0.25*t) - 0.25*t;
    Hst[ 2] = -0.125*r*r - 0.25*s + r*(-0.125 - 0.25*s + 0.25*t) + 0.25*t;
    Hst[ 3] = -0.125*r*r - 0.25*s + r*(0.125 + 0.25*s - 0.25*t) + 0.25*t;
    Hst[ 4] = 0;
    Hst[ 5] = -0.25*(-1. + r*r);
    Hst[ 6] = 0.5*(1. + r)*s;
    Hst[ 7] = 0.25*(-1. + r*r);
    Hst[ 8] = -0.5*(-1. + r)*s;
    Hst[ 9] = -0.5*(-1. + r)*t;
    Hst[10] = 0.5*(1. + r)*t;
    Hst[11] = -0.5*(1. + r)*t;
    Hst[12] = 0.5*(-1. + r)*t;
    
    Hrt[ 0] = r*(-0.25 + 0.25*s) + 0.125*s*s + s*(-0.125 + 0.25*t) - 0.25*t;
    Hrt[ 1] = r*(-0.25 + 0.25*s) - 0.125*s*s + s*(0.125 - 0.25*t) + 0.25*t;
    Hrt[ 2] = r*(-0.25 - 0.25*s) - 0.125*s*s + s*(-0.125 + 0.25*t) + 0.25*t;
    Hrt[ 3] = r*(-0.25 - 0.25*s) + 0.125*s*s + s*(0.125 - 0.25*t) - 0.25*t;
    Hrt[ 4] = 0;
    Hrt[ 5] = -0.5*r*(-1. + s);
    Hrt[ 6] = 0.25*(-1. + s*s);
    Hrt[ 7] = 0.5*r*(1. + s);
    Hrt[ 8] = -0.25*(-1. + s*s);
    Hrt[ 9] = -0.5*(-1. + s)*t;
    Hrt[10] = 0.5*(-1. + s)*t;
    Hrt[11] = -0.5*(1. + s)*t;
    Hrt[12] = 0.5*(1. + s)*t;

    Hrs[0] = 0.125*(1.0 - t); Hrt[0] = 0.125*(1.0 - s); Hst[0] = 0.125*(1.0 - r);
    Hrs[1] = -0.125*(1.0 - t); Hrt[1] = -0.125*(1.0 - s); Hst[1] = 0.125*(1.0 + r);
    Hrs[2] = 0.125*(1.0 - t); Hrt[2] = -0.125*(1.0 + s); Hst[2] = -0.125*(1.0 + r);
    Hrs[3] = -0.125*(1.0 - t); Hrt[3] = 0.125*(1.0 + s); Hst[3] = -0.125*(1.0 - r);
    Hrs[4] = 0.0; Hrt[4] = 0.0; Hst[4] = 0.0;
}
